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Abstract 

In this paper, we study degenerate Hopf bifurcations in a class of parametrized 
retarded functional differential equations. Specifically, we are interested in the case 
where the eigenvalue crossing condition of the classical Hopf bifurcation theorem is 
violated. Our approach is based on center manifold reduction and Poincare-Birkhoff 
normal forms, and a singularity theoretical classification of this degenerate Hopf bi- 


1 



furcation. Our results are applied to a recently developed SIS model incorporating a 
delayed behavioral response. We show that the phenomenon of endemic bubbles, which 
is characterized by a branch of periodic solutions which bifurcates from the endemic 
equilibrium at some value of the basic reproduction number Rq, and then reconnects to 
the endemic equilibrium at a larger value of Rq, originates in a codimension-two orga¬ 
nizing center where the eigenvalue crossing condition for the Hopf bifurcation theorem 
is violated. 
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1 Introduction 


Retarded functional differential equations (RFDEs), of which delay differential equations 
(DDEs) are a special case, are used to model a large variety of phenomena in sciences, 
engineering, economics and many other areas [H El |3l IH EHl 1121 |l3l ESI ESI ESI ESl EDI El] • 
One of the main technical differences between RFDEs and ordinary differential equations 
(ODEs) is that while ODEs require initial data only at one point (typically at time t = 0) 
to generate a solution, RFDEs require initial data in a range of past values, typically an 
interval t G [—r, 0]. One of the consequences of this fact is that the phase space for RFDEs 
is inhnite-dimensionah Despite this fact, most of the usual tools and techniques of the theory 
of dynamical systems can be suitably adapted to the study of RFDEs [H E] IH|. In particular, 
for parametrized RFDEs, one can analyze bifnrcations using center manifold theory and 
normal forms BUS], and in many studies (see for example [21 El EDI ESI ESI ED]), this has led 
to valuable insight into many phenomena which are modeld using RFDEs. 

The focus of this paper is an analysis of a certain codimension 2 degenerate Hopf bifnr- 
cation in RFDEs. In particular, we are interested in stndying the case where a parametrized 
RFDE admits an eqnilibrium solntion which satishes the simple purely imaginary eigenvalue 
condition of the Hopf bifnrcation theorem |13], but for which the crossing condition of this 
theorem is violated. This is motivated from recent results m where a SIS model incor¬ 
porating a delayed behavioral response was analyzed and shown to exhibit a phenonmenon 
which the authors called endemic bubbles. It became apparent to us that the bifurcation 
diagrams reported in m resembled bifnrcation diagrams which had been fonnd to be in 
the versal unfolding of the degenerate Hopf bifurcation (with crossing condition violated) in 
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[B]. Although this SIS model is the application which motivated our analysis, the theoretical 
results we present here have a much larger scope of application, since many phenomena in 
nature are modelled using the class of delay differential equations we study. 


1.1 The crossing condition for the Hopf bifurcation 

Consider as a prototype the delay differential equation 


x{t) = ax{t) + {3 x(t — r) + F{x(t),x(t — r)) 


( 1 . 1 ) 


where a, f3 E M. are parameters, r > 0 is a hxed delay time, and F represents higher-order 
nonlinear terms. 


Equation (1.1) has x = 0 as an equilibrium solution for all values of a and /9. We are 


interested in bifurcations from this trivial equilibrium point. To this end, we consider the 
characteristic equation 

^ = a + 13 

In particular, we are interested in purely imaginary solutions ^ = iu io this characteristic 
equation. It is easy to see that such solutions occur when the parameters a and [3 satisfy 


a + (3 cos TOJ = 0, —uj = 13 sin rcn, uj = \/(3“^ — (3"^ > c?. 


( 1 . 2 ) 


Equations (1.2) define curves in the a-/? parameter space, as illustrated in Figure 


Now let us suppose that a and [3 in (1.1) depend on a distinguished external control 


parameter A. If the path (q!(A),/ 9(A)) in parameter space crosses the Hopf bifurcation curve 


(1.2) at the point [a*, (3*) = (a(A*),/9(A*)) from the region where the equilibrium x = 
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Figure 1; Hopf bifurcation curves (1.2) in red for r = 1 on the left and r = 5 on the right. 
The blue line represents (3 = a. In the region in parameter space below the blue line and 


above the uppermost red curve, the equilibrium point x = 0 of (1.1) is locally asymptotically 
stable, and loses this stability below the uppermost red curve. 
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0 is stable into the region where it is unstable as illustrated in Figure then typically 
it crosses this Hopf curve transversally, and assuming non-degeneracy conditions in the 
nonlinear coefficients (specihcally, the hrst Lyapunov coefficient is non-zero), a super- or 
sub-critical Hopf bifurcation occurs from the trivial equilibrium. 



Figure 2: As the control parameter A varies through A*, the (black) curve (q!(A),/9(A)) 


crosses the (red) Hopf bifurcation curve (1.2) at the point (a*,/9*). Generically, the crossing 
is transversal (as illustrated here), and we get a Hopf bifurcation from the trivial equilibrium 


point of (1.1). 


In applications, it often occurs that in addition to depending on the distinguished pa¬ 
rameter A, the coefficients a and f3 in (1.1) also depend on auxiliary parameters p G MF. For 
the purposes of this discussion, let us suppose that p = 1. In this case, it is possible that for 
a certain value of /i = p*, the curve A i —> (a(A, p*), p*)) intersects the Hopf curve (1.2) 

tangentially at the point {a*, (3*) = {a{X*, p*), (3{X*, p*)) as is illustrated in Figure]^ In this 


6 







case, the transversality condition of the Hopf Bifurcation Theorem [T3] is violated. However, 
as fi varies in a neighborhood of /r*, this degeneracy is “unfolded” as illustrated in Figure 
1^ Assuming certain non-degeneracy conditions on nonlinear terms for ODEs undergoing a 
violation of the crossing condition, it is shown in ^ that the Hopf bifurcation diagrams in 
a neighborhood of fi* are equivalent to one of the bifurcation diagrams of the normal form 

x{e{X‘^ + 7]) + x‘^) = 0, e = ±1 (1.3) 



Figure 3: Tangential intersection of the parameter curve (a(A,/r*),/3(A,/i*)) with the Hopf 


bifurcation curve (1.2) at A = A*. The crossing condition of the Hopf bifurcation theorem is 


violated for this parameter path. 


Consider the twice differentiable path A i —> (a(A,/r),/9(A,/r)), and let A* and p* be such 
that the point 

(«*,/?*) ^(«(A*,/x*),/3(A*,/i*)) 
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Figure 4: Unfolding the tangential intersection (middle diagram), for fixed values of /i in a 
neighborhood of /i*. Generically, we either get no intersections (as illustrated on the left) or 
two intersections (as illustrated on the right) with the Hopf bifurcation curve. 


satishes both the following equations 

a* + f3* cos r a/= 0 


(1.4) 

(1.5) 


then the path A i —> (a(A, fi), /3{X, /i)) has a tangential intersection with the Hopf bifurcation 


curve (1.2) at the point (a*,/?*). At the point (a*,/3*), the signed curvature of the Hopf 


bifurcation curve (1.2) can be computed as 

- /?*2)(/5*V2 + a*T - 2)r 


Kl = 


[(/S*^r^ + l)(a*^ + ( 3 *^) — 4a*/S*^r]^A 


and the signed curvature of the path A 

{a*r - - 

/3x{X*, 


I—^ (a(A,/i),/3(A,/i)) is 

- l)axxiX\i^*) + (a* - P*^t)I3xx{X\ii*)) 
+ l)(a*^ + /d*2) — 4a*/d*2r]^A 


K2 = 











Figure 5: Universal unfolding of the degenerate crossing condition for the Hopf bifurcation. 


These are the bifurcation diagrams for the normal form (1.3). An “s” designates a sta¬ 
ble branch of periodic solutions, whereas a “u” designates an unstable branch of periodic 
solutions. 
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Although these two curves are tangential at the intersection point, we want their curvatures 
to be different, so that locally near [a* , (3*) the curve A i —> (a(A,/r),/9(A, p)) lies entirely on 
one side only of the Hopf curve, as illustrated in Figure This is equivalent to requiring 

g(a*, /?*, r) = {a*T - if ^J*) + - if {a* - /3*V)/3aa(A*, p*) 

( 1 . 6 ) 

-(/?a(A*, /i*))V(a*2 _ + a*T -2) fO. 

By using center manifold and normal form analysis for retarded functional differential 
equations (see for example 010). we expect to see for (A,/x) near (A*,/i*) Hopf bifurcation 
diagrams for 

x(t) = a(A, /i) x(t) + (3{\, /i) x(t — r) + F{x(t), x(t — r). A, /i) 
which resemble those of Figure]^ This will be the theoretical focus of section 2 of this paper. 

1.2 Endemic bubbles 

In m, the following susceptible-infected-susceptible (SIS) model incorporating delayed be¬ 
havioral response was analyzed 

yf ) = -y{t) + Roh{y{t - T))y{t){l - y{t)) (1.7) 

where y{t) represents the proportion of infected individuals in the population, Rq is the basic 
reproduction number (expressing the expected number of secondary infections generated by 
a single infectious agent introduced into a wholly susceptible population), and the smooth 
behavioral function h : [0,1] — )■ (0,1] is such that h'{y) < 0, h(0) = 1 and h(l) < 1. The 
parameter r > 0 represents the delay in time between the moment when the population 
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has information about the number of infections, and the moment this population decides to 
modify its behavior in order to reduce the rate of infections. We refer the reader to m for 
more details related to the modelling aspects. 


Equilibria of (1.7) are algebraic solutions to 


y = Pu3h(y)y(l-y). 


( 1 . 8 ) 


Obviously ?/ = 0 is always a solution to (1.8) (disease-free equilibrium). It is shown in 
m that if i?o < then the disease-free equilibrium is globally asymptotically stable. If 


Rq > 1, then y = 0 becomes unstable, and (1.8) has a unique endemic equilibrium y 


satisfying y* < 1 — 1 /Rq. The authors in [TT] then perform a comprehensive analysis of 


(1.7) for different choices of response functions h, and in particular, they plot bifurcation 


diagrams with distinguished parameter Ro. Some of these bifurcation diagrams (see for 
example hgures 5 and 9 of m) exhibit a phenomenon that the authors have called endemic 
bubbles, which loosely speaking is the bifurcation of a branch of periodic solutions from the 
endemic equilibrium point at some parameter value i?o > 1, and this branch reconnects with 
the endemic equilibrium (in a reverse Hopf bifurcation) at parameter value Rq > Rq. A 
schematic representation of a typical such endemic bubble bifurcation diagram is given in 
Figure 

When focussing on the endemic branch, we notice the similarity between Figure and 
the bottom-right-most bifurcation in Figure]^ As we will show in this paper, the endemic 


bubble is indeed a consequence of a degenerate Hopf bifurcation which occurs in (1.7) for 


various choices of response functions h{y,p), where p > 0 is some auxiliary parameter which 
appears in the behavioral response function. 
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Figure 6: Schematic representation of a bifurcation diagram for (1.7) which exhibits the 
phenomenon of endemic bubble. The disease-free equilibrium is stable for Ro < 1. At Rq = 1, 
a branch of endemic equilibria bifurcate from the disease-free equilibrium, and this branch 
is stable until Ro = Ro at which point a Hopf bifurcation occurs. The blue curves represent 
minima and maxima of oscillations for the bifurcating periodic solutions. At Ro = Ro, the 
branch of periodic solutions reconnects with the endemic equilibrium. The region between 
Ro and Ro is what is referred to as the endemic bubble. 
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1.3 Overview of paper 


With the above discussion in mind, this suggests that the organizing center of endemic bub¬ 


bles is a degenerate Hopf bifurcation in (1.7). In this paper, we will develop the theoretical 


ideas to make this hypothesis into a rigorous theorem. Our approach will be to use the 
center manifold and normal form theory for RFDEs developed in [U |5] and the classihcation 
and unfolding of degenerate Hopf bifurcations of [B] . 

Ours is not the hrst study of higher codimension bifurcation in RFDEs with nonlinear 
degeneracies. We note in particular [2] where the authors present and analyze a DDE model 
for tumor growth in which a Bautin bifurcation (Hopf bifurcation with the hrst Lyapunov 
coefficient being zero) occurs. In na, a Maple program is presented which allows the com¬ 
putation of coefficients up to any order related to Hopf bifurcation for RFDEs. Also, in [3] , 
the authors investigate the stabilization of unstable periodic orbits via the Pyragas delayed 
feedback control. It is shown that the stabilization mechanism occurs because of a highly 
degenerate Hopf bifurcation in which both crossing condition and hrst Lyapunov coefficient 
degeneracy occur simultaneously. Our paper here appears to be the hrst systematic study of 
the degeneracy resulting from the violation of the crossing condition in the Hopf Bifurcation 
Theorem for general RFDEs. 

In section 2, we will establish the sufficient conditions for such a degenerate Hopf bifur¬ 
cation to occur in a general class of parametrized RFDEs. Our main theorems to that ehect 


are Theorems 2.3 and 2.4 In subsection 2.3, we will apply these results to the special case 


in which the RFDE is a DDE of the form (1.1) 


In section 3, we use our theoretical results to seek out and analyze degenerate Hopf 
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bifurcation points in the SIS model (1.7), where the behavioral response function h = h{y,p) 


depends on an auxiliary parameter p. The parameter space is thus two-dimensional: Rq acts 
as a distinguished bifurcation parameter, and p as an unfolding parameter. This leads to 


Theorem 3.1 In subsections 3.2 and 3.3, we illustrate these results by performing explicit 
computations for two of the types of delayed response functions h{y, p) which were considered 

in [TT] . 

We end with some concluding remarks in section 4, and an Appendix where a lengthy 
expression for the Lyapunov coefficient associated with the degenerate Hopf bifurcation is 
presented. 


2 Center manifold and normal form reduction 

In this section, we give a brief summary of the theory presented in mu for center-manifold 
and normal form reductions of retarded functional differential equations, and then apply it 
to study the degenerate Hopf bifurcation in a class of parametrized RFDEs. 

2.1 Phase space and splitting 

For r > 0, we consider the Banach space C = C([—r, 0],M) of continuous functions from 
[—r, 0] into M, endowed with the supremum norm. We dehne Zt & C hj Zt{9) = z{t + 9), for 
9 e [-r,0]. 

Let 

z{t) = L{X,p){zt) + F{zt,\,p) (2.1) 


14 




denote a smoothly parametrized family of nonlinear retarded fnnctional differential eqna- 
tions, where A G M is a real distinguished bifurcation parameter, /r G M is an auxiliary 
parameter (to be regarded as an unfolding parameter), L(A,/i) is a smoothly parametrized 
family of bounded linear operators from C into M, and F is a smooth function from C x 
into M, which satisfies 


F(0, A,/i)=0, iAiF(0, A, p) = 0, V(A,/x) near (0,0) G 


( 2 . 2 ) 


We denote Lq = L(0,0), and we rewrite (2.1) as 


z{t) — Lo^Zt) + {L{X, fi) — Lo)(xi) + F{zt, A, /i). 


By the Riesz representation theorem, we may write 

= j Na,m(^)]<(>(^), 

where ri\^^{6) is a measurable function on [—r, 0]. We also define Al(A, /r) to be the infinites¬ 
imal generator for the flow of the linear system z = L{X, ^){zt). For each fixed (A,p), we 
consider the characteristic equation 

A(A, ^)(0 = 0, A(A, ^)(0 dmje) e«». 

We suppose that the following holds: 

Hypothesis 2.1 The characteristic equation A(A,/r)(.^) = 0 has roots, ^(A,/i) = 7 (A,/i) ± 
icj{X,fi) which are smoothly parametrized by X and /i, and are such that 

7(0,0) = 0, a;(0,0)=a;*>0, 7 a(0,0) = 0, 7/a(0,0)^0, 
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(2.3) 


—A(0, 0)(za;*) 7 ^ 0, or equivalently 1 — Lq{6 7 ^ 0, 


and furthermore, we suppose the characteristic equation A(0,0)(^) = 0 has no roots other 
than ±iu* on the imaginary axis, and that all other roots of = 0 have strictly 

negative real part. 

It follows that if Aq denotes the inhnitesimal generator of the linear equation z{t) = 
Lq^Zi), then the eigenspace P of yio corresponding to the eigenvalues dnioj* is two-dimensional. 
Let *h( 6 ') = be a complex representation of a basis for P. The dual space 

C* = C([ 0 , r], M) also admits a two-dimensional subspace for the transpose Aq corresponding 
to dziuj*. We introduce the bilinear form between C and C* 



We choose a basis \h(s) = col(V’i(0)e '0i(O)e*‘^®) such that (th, <h) = Id- As is shown in 
[5] we have 




(2.4) 


which is well-defined because of (2.3). We have C = P ® Q, where Q is infinite-dimensional 


and also invariant for the operator ^Iq. 

We now let BC represent the Banach space of functions from [—r, 0] into M which are 
uniformly continuous on [—r, 0) with a jump discontinuity at 0. We define the function 




1 , 0 = 0 
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The elements of the space BC can be written as F = (^ + where ip G C and z/ G M. We 
dehne the projection tt : BC —)■ P as 

7r{(p + Xqu) = <h[(T, if) + \I>(0)z/]. 


We then have 


BC = P® her vr, 


where Q C kervr, and we may rewrite (2.1) according to this splitting as 


X = Bx + \['(0)[(L(A,/r) - Lo){^ x + () + F{^x + C, A,p)], 

|C = ^o.C + (h - ir)W[(i(A, A.) - LM<S> x + C) + F(<S>x + <;, A, ft)], 


(2.5) 


where 


B = 


^ iu* 0 ^ 


V 


/ 


0 -iu* 

X = {u, ^ = Q n C^([—r, 0], M), and Aq^ is defined as 


Aq^ip = ip + Xo[Loip - 0 ( 0 )]. 


As is shown in |5], it is then possible to define near identity changes of coordinates 


(a:,C) = (£,C) + U{x,\p) 


( 2 . 6 ) 


such that a Taylor expansion of (2.5) can be put into the normal form 


d 


X = 


Bx + Y^ ^)(a;,C, A,/x), = ^QiC + Y 


i>2 


i>2 
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having the property that the center manifold is locally given by ^ = 0, and the local flow of 


(2.1) on this center manifold is given by 


X = 


Bx + Y^ g]{x,0,X,n). 

i>2 


(2.7) 


The nonlinear terms in (2.7) are in classical Poincare-Birkhoff normal form with respect to 
the matrix B. 


2.2 Normal form for the degenerate Hopf bifurcation 


From Hypothesis 2.1, we have that there exists a smooth function ,^(A,/i) such that 


e(A,/i) = L(A,/i)(e«("’'^)'), e(0,0)=za;*, 


for which implicit differentiation gives the following equalities 


a(0,0) =^/>i(0)LA(0,0)(e 






and 


5aa(0,0) = i/.,(0) [LAA(0,0)(e”’'') + 2J,(0,0)LA(0,0)(9e''"'') + KA(0,0))V„(9'^e-“‘")] , 

( 2 . 8 ) 


where V’i(O) is as in (2.4). Hypothesis 2.1 thus implies 


Re [V'i(0)LA(0,0)(e*^*'')] = 0, and Re [V’i(O) ^ 0. 


In addition, we are going to assume the following second-order non-degeneracy condition 
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Hypothesis 2.2 We assume the non-degeneracy condition Re(^AA(0,0)) 7 ^ 0, which is equiv¬ 


alent (via (2.8)) to 


Re{M0) [LAA(O,O)(e*‘^*^) + 2a(O,O)LA(O,O)(0e*‘"*^) + (eA(O,O))%(02gi.*e)]) ^ ^2.9) 


We will see later on that the geometrical meaning of (2.9) is equivalent to the curvature 


condition given in (1.6). 


Consider equations (2.5) which we expand in a Taylor series as 


X = 5x +T(0 )[ALa(0,0)($x + C)+/i-h;,(0,0)(<l)a: + C) 


dt 


+ 2'^^ Laa( 0, 0)($a: + C) + 


C — ^QiC + (-^rf ~ -^a( 0,0)($X + C) +/^-h^(0,0)(<hx + 0 


( 2 . 10 ) 


+ 2-^ -^aa( 0 , 0 )(<hX + C) +-R(a^,C)-^5h )]5 
where the remainder term is /2(x,C, A,/r) = 0(|(x, /i^, pA, A^). 

We precede as in section 3 of [5]. Writing x = {u,u)'^, $x = we dehne 


quadratic near-identity changes of coordinates of the form (2.6) which transforms (2.10) in 


such a way that the normal form (2.7) becomes 


ii = iu*u + [((Ti -h i(T 2 )ia + ia^X ( 0-4 ia^)X^]u + V’i(0)i?((M, m), 0, A,p), 


( 2 . 11 ) 


where 

(Ti ia2 = V’i(O) T^(0, 0)(e*‘^*®), ai ^ 0, 

as = Im [V’i(0)TA(0,0)(e'‘"‘^)] = a;A(0,0), 
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and 


a4 + *a5 = IMO) [i:;,,(O,O)(e-*^) + 2eA(O,O)LA(O,O)(0e-*^) + (eA(O,O))2Lo(02e-‘®)] 


= kAA(0,0) 


( 2 . 12 ) 


is such that <74 7 ^ 0 because of Hypothesis |2.2[ The extra terms 


-V^i(O) [ 2 eA(O,O)LA(O,O)( 0 e“"*^) + {^0, 0)^10(6^ e^‘ 


iuj*6 


which appear in ( 2 . 12 ) arise from having normalized quadratic terms in ( 2 . 10 ). 


As mentioned earlier, the nonlinear terms in (2.7) (equivalently ( 2 . 11 )) are in classical 
PoincarABirkhoff normal form with respect to the matrix B. It is well-known (see for 
example 0 ) that such a normal form has the algebraic form 


u = iu*u + r(|M|^, A, jj)u^ 


(2.13) 


where P is a smooth complex-valued function. 

The following theorem now follows from the above discussion and by performing further 


near-identity changes of coordinates (2.6) as in [5]: 


Theorem 2.3 Consider the smoothly parametrized family of nonlinear retarded functional 


differential equations (2.1) which satisfies (2.2), and Hypotheses 2.1 and 2.2 Then there 
exists a two-dimensional semi-flow invariant center manifold in the phase space C'([—r, 0], M). 
Furthermore, there exist a formal sequence of parameter-dependent near-identity changes of 


coordinates of the form (2.6) which are such that the dynamics of (2.1) reduced to this center 


manifold have a Poincare-Birkhoff normal form to any order given by (2.1) (equivalently 
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(2.13)): 


u= [icu* + {ai + ia2)fi + ia3X + {a4 + ia5)X‘^ + H{\u\‘^,X,fi)]u, (2.14) 

where H{\u\'^, 0 , 0) = Ki + iK 2 , ifAd^P, 0 , 0 ) = 0 , 0 ) = Hxx{\u\‘^, 0 , 0 ) = 0 . 


The coefficient Ki in (2.14) is the hrst Lyapunov coefficient. It has been computed explicitly 
in [S] as 


Ki = Re 


„/, rn\ I D ^(1,1,0,0)5(1,0,1,0) , 5(2,o,0,0)5(o,1,0,1) 

^l(U) 13(2,1,0,0)--1- 


^o(l) 


2iuj* — 


(2.15) 


where the coefficients are read-off from the Taylor expansion of F in (2.1): 

F(a;ie*‘^*® - 1 - X2e~^‘^*^ -t- X 3 1 - 1 - 0 , 0 ) = 

5(2,0,0,0)3^1 + 5(i^i_o,0)3;ia:2 + 5(i_o,l,0)^1^3 + 5(oy,o,l)^2a^4 + 5(2,1,0,0)3^13^2 + ■ ■ ■ 


(2.16) 


Using the polar coordinates u = (2.14) becomes 


r = r((j 4 A^-I-(Ti/i-I-iLir^-I-0(r^, pA, A^)) 


(2.17) 


= u* - 


Zeros of the r equation of (2.17) correspond to periodic solutions of (2.1) whose period is 


approximately equal to 2ti/u* . Let us assume the non-degeneracy condition TLi 7 ^ 0. Since 
cTi 7 ^ 0 and CT4 7 ^ 0, it now follows from the classihcation of degenerate Hopf bifurcations 
done in [ 6 ] that 

Theorem 2.4 For fi near 0, the bifurcation diagram (wrt X) of r(a 4 X‘^ + cxi/i -|- Kir^ + 
0 (r^p^/iA, A^)) = 0 is locally diffeomorphic to the bifurcation diagram of 


r{e{X‘^ + h) + = 0) 
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where e = sgn((T 4 /i^i) and rj = 


CTi/i 


|/Ci|sgn(cr4) 


(see (1.3)), as illustrated in Figure 


5 


2.3 Special case: discrete delay 

As a special case of the theory we have just developed, we return to the prototype delay- 


differential equation ( 1 . 1 ) 


x{t) = q;(A, ti)x{t) -|- /3(A, fi)x{t — r) -|- F{x{t),x{t — r), A, /i) 


(2.18) 


where we will assume without loss of generality that A* = 0, p* = 0, and {a*, (3*) are such 


that (1.4), (1.5) and (1.6) are satished, and (3* ^ 0, and — 1 7 ^ 0. 


The operator L{\, fi) is dehned as 


L(A,/i)(z(0)) = a{\, fr)z{0) + f3{\, ix)z{-t), 


so that 


Lo(e*‘"^) = a* + (3*e 


^ = iuj* ^ e-^^ ^ 


(3* 


Lo( 0 e*‘" = -T(3*e-^^ " = -riiuj* - a*), 


^i(O) = [i-Lo{eF^*%]-^ = 


Using (1.5), we have 

V^i(0)L,(0,0)(e-*^) 




(1 — a*T) + ioj*T 

«A(0,0) + ^^(*n;*-«*) 

(1 — a*r) -f iuj*T 

n;*(«A(0,0)/3*r-/3A(0,0)) 


/d*((l — Q;*r)^ -|- 

The non-degeneracy condition Re('0i(O)L^(O, 0)(e*‘^*®)) becomes 

_ /3*«^(0,0)(1 - a*T) + /3^(0, 0){r/3*^ - a*) 


CTi = 


/d*((l — Q;*r )2 -|- 


^ 0 . 
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Finally, 

2(74 = Re(V^i(0)[L,,(0,0)(e-*^) + 2^7(0, O)LA(O,O)(0e-*^) + (^ 7 ( 0 , O))2Lo(02 

_ 5{a*,/3*,T) _ 

f3*'^{a*T — + 2a*T — 1) ’ 


where ^{a* ^ j3* ^r) is as in (1.6). Thus, we note that the condition <74 7 ^ 0 is equivalent to 


the curvature condition ( 1 . 6 ). 


We write F in (2.18) at (A,/i) = (0,0) as 
F{x{t),x{t - r), 0, 0) = + f(i^i)x{t)x{t - r) + f(o, 2 ){x{t - T)f+ 

+f{3,o)x{tf + f(^2,i)x{tfx{t - r) + /(i,2)x(t)(a;(t - rjf + f{o,3){x{t - r))' 


+0(|a;| 


(2.19) 


To compute the first Lyapunov coefficient Ki in (2.15), we need 


{ico* - a*f 
f3* 


Lo(l) = a* + /3*, and Lo{e^^^ ^) = a* + 13*e 


D* ^—2iu*r 


* I 

— a T 


We will assume that a* +^* 7 ^ 0. If this condition holds, it follows that 2ioo* — Lo{e^''‘^ ®) 7 ^ 0. 


A lengthy computation using (2.15), (2.16) and (2.19) yields the Lyapunov coefficient Ki in 


terms of the coefficients /(y/c) in (2.19) and a*, /?*, uj* and r. The formula is lengthy and 
given in the Appendix. Generically, we will have iFi 7 ^ 0. Once we have these quantities, 
we can compute the normal form for the unfolding of the degenerate Hopf bifurcation at 
(A,/i) = (0,0): 

r(e:(A^ + h) + = 0) 


where £ = sgn((T 4 /iFi) and rj = 


CTi/X 


|Wi|sgn((T4) 


as in Theorem 


2.4 
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3 Application to a SIS model with delayed behavioral 


response: endemic bubbles 


We recall the SIS model (1.7) which was analyzed in m- In this section, we will apply the 
theoretical resnlts of the previous section to establish and study degenerate Hopf bifurcations 
which occur in this model. 


3.1 General case 


Non-trivial equilibria of (1.7) are solutions y 0 to (1.8), or equivalently to 

h{y,p) = 


(3.1) 


fio(i - y)' 

Because the function h is such that hy{y,p) < 0, h{0,p) = 1 and h{l,p) < 1, if i?o > 1 then 


(3.1) has a unique solution y = y{Ro,p), i.e. 


Roh{y{Ro,p),p){l - y{Ro,p)) = 1, Vi?o > 1,P > 0 . 


Linearizing (1.7) about the equilibrium y{RQ,p) gives 


where 


and 


x{t) = a{Ro,p)x{t) + (3{Ro,p)x{t - r). 


a{Ro,p) = -Roy{Ro,p)h{y{Ro,p),p) < 0 


(3.2) 


(3.3) 


(3{Ro,p) = Rohy{y{Ro,p),p)y{Ro,p) (1 - y{Ro,p)) < 0. 


(3.4) 
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Using implicit differentiation of (3.2), one can compute yp and in terms of h, hy, 


hyy, hp and y. One can then use (3.3) and (3.4) to compute the quantities /3p, 

C^RoRo PrqRq. 

We then have the following 


Theorem 3.1 Consider the nonlinear delay-differential equation 


x{t) = a{Ro,p)x{t) + (3{RQ,p)x{t - r) + 3^{x{t),x{t - r), Ro,p) 


which is obtained by performing the change of variables x = y — y{Ro,p) in (1.7), (where 


y{Ro,p) is the endemic equilibrium, which is solution to (3.2)), and where a{Ro,p) and 


/3{Ro,p) are as in (3.3) and (3.4) respectively. Define u}{Rq,p) by 


uj{Ro,p)‘^ = /3(i?o,p)^ - a{Ro,p)‘^. 


Suppose {Rq,p*) is a point in parameter space such that 


a*(3* cos u*T = 0, 


(3.5) 




(3.6) 


where a* = a{Rl,p*), fi* = fi{R*Q,p*), a*^ = aRfiRl,p*), (3*^ = fi{Rl,p*), a* = ap(Ro,p*), 
(3; = (3p{Rl,p*), a)sgKg = aR,RfiR*o,p*), (3),^^^ = fiR,RfiR*„p*), and = (3*^ - > 0. 

Suppose 

(3*0* il - a*T'\ + B* (tB*"^ - oO 

^ 0, (3.7) 


, ^ /3*a;{l-a*T)+fi;{T/3*^ -a *) 
/3*{{1 — a*T)‘^ 

2 _ 9{a*,(3*,r) _ 

^ (3*‘^{a*T — l)‘^{—f3*‘^T‘^2a*T — 1) 


7 ^ 0 , 


(3.8) 
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(where S is as in (1.6)), and that the first Lyapunov coefficient K( = Kfia*, (5*) computed in 


Appendix A is non-zero. Then 0 has a degenerate Hopf bifurcation point at the endemic 
equilibrium y when {Ro,p) = {Rq^P*)) normal form for the Hopf bifurcation diagrams 

near this point is 

r{e{\^ + 7]) + r^) = 0 

p — p*) 

\sgn{al 

enough to {Rq,p*), this Hopf bifurcation diagram will exhibit an endemic bubble if e = -\-l 
and rj < 0. In the case e = +1 and rj < 0, the width of the endemic bubble (in Rq space) is 
approximately equal to the width of the region between the two zeros of 

afiRo - Rq)"^ + afip - p*), 


Furthermore, for {Ro,p) close 


where X = Rq — Rq, e = sgffiaX/Kfi) and rj = 


\Kr 


i.e. 


Ra — Rn 


afip-p*) 


(74 


3.2 Response function h{y,p) = 

This is one of the special cases which was studied in m and for which endemic bubbles were 
observed. For purposes of comparing our results to those of m, we will assume as they do 
that T = 10. 

The endemic equilibrium is 

Rq — 1 


y{Ro,p) = 


p + R 


i?o > 1 
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and introducing x = y — y* transforms (1.7) into 


x{t) = a{Ro,p)x{t) + (3{Ro,p)x{t - r) + 0{\x\'^), 


(3.9) 


where 


/D X 1-^0 ..p ^ (l-^o)p 

a{Ro,p) = — -, /3{Ro,p) = 


1+p' Ro{l+p)' 

Computing the quantity — ar) + Pro{t/ 3^ — a) for r = 10 gives 

(1- Rofp (lO Ro^ + Ro^p + Ro^ - 10 p^) 
Ro\i + pf ■ 


So using (3.6) and solving for p, we have 


p — p* — — ( i?o + V -^0 + 400 i?o + 40 ) i?0) 


where Rq is a root of (3.5), which is equivalent to 


where 


= 2 


cos (10 w*) + 400i?S + 40^ + 20 = 0 


(1 - Rif (rI ^/Rf + 400 Rl + 40 + Rf + 200 Rl - 180 


(3.10) 


Rl ^/Rf + 400 Rl + 40 + Rf + 20 
A plot of the left-hand side of (3.10) as a function of Rq is given in Figure]^ We numerically 
compute the value for the root 


Rl ^ 1.784 


which yields 


p* ^ 2.613 
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Figure 7: Graph of the left-hand side of (3.10). 


and 


a 


-0.217, ^ -0.318, (u* ^ 0.232. 


We compute and as in (3.7) and (3.8) and get 


at ^ 0.021, at ^ -0.037. 


To compute the Lyapunov coefficient iFf, we will need to compute explicitly the quadratic 


and cubic terms in (3.9). A straightforward computation gives the following expression for 








these quadratic and cubic terms: 




(p*+-R5)(i?S-l)p* 


Rfil+P* 




+ (%Tiw) - ^) - P*) ) x{t)x{t 

The formula for Ki in the Appendix now gives 




so that 


Kl ^ -1.006 


£ = sg\i{aX/Kl) = +1, 1] ^ —0.021 (p — p*). 


Based on Theorem 3.1, we therefore predict the following for p near p* when h(y,p) = 
1/(1+ yp): 


(i) Suppose p < p*. Then varying Rq near in (1.7) will not generate Hopf bifurcation 
from the endemic equilibrium, hence there will be no endemic bubble. 


(ii) Suppose p > p*. Then varying Rq near R^ in (1.7) will generate an endemic bubble 
via Hopf bifurcation at Rq = Rq followed by a reverse Hopf bifurcation at Rq = Rq, 
with Rq and Rq near Rq, and 


-Rq ~ -Rq 


a 


— (p — p*) ~ 1.614i/p — p*. 


(Ta 


In Figurewe illustrate the results of several integrations of (1.7) with h{y,p) = 1/(1 + 
py) for values of {Rq, p) near {Rq, p*). The resulting bifurcation diagrams show the emergence 
of the endemic bubble. 
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{a) p=2.61 


(b) p=2.62 


y 0.1 


1.5 1.6 1.7 1.8 1.9 2 2.1 

(c)p=2.633 


***** 

* Jp 


y 0.18- 


□ O 

O 

o qO 

^ OoO° 


1.5 1.6 1.7 1.8 1.9 2 2.1 


•k r-,n 


□° OoO° 


1.5 1.6 1.7 1.8 1.9 2 2.1 


(d)p=2.7 


^******* 


0 . 26 - 

0 . 24 - 

0 . 22 - 

0 . 20 - 

0 . 18 - 

0 . 16 - 

0 . 14 - 




o 

0.12ti]° o „0 

°Onor>nO°° 


1.5 1.6 1.7 1.8 1.9 2 2.1 


Rr, 


Figure 8: Numerically obtained bifurcation diagrams for (1.7) with h{y,p) = 1/(1 + py) 


for the following values of p: (a) p = 2.61 < p*, (h) p = 2.62 > p*, (c) p = 2.633 > p*, 
{d) p = 2.7 > p*. The red squares indicate that the simulation has settled to the endemic 
equilibrium. The green asterisks and blue circles designate respectively the maxima and 
minima of the steady oscillations which occur for values of Rq inside the endemic bubble. 
Recall that the theoretically computed value of {Rq,p*) for the degenerate Hopf bifurcation 
is {Ro,p*) ^ (1.784,2.613). 
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(a) p=1.660 


y 0.27 


2.1 2.2 


(b)p= 1.662 


y 0.27 




Figure 9: Numerically obtained bifurcation diagrams for (1.7) with h{y,p) = e ™ for the 


following values of p: (a) p = 1.660 < p*, (b) p = 1.662 > p*. The red squares indicate 
that the simulation has settled to the endemic equilibrium. The green asterisks and blue 
circles designate respectively the maxima and minima of the steady oscillations which occur 
for values of Rq inside the endemic bubble. Recall that the theoretically computed value of 
{Rq,P*) for the degenerate Hopf bifurcation is {Rq,p*) (2.1474,1.6617). 
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Figure 10: Numerically obtained periodic solution y(t) to (1.7) with h{y,p) = e~^ after 
transients have died out. The parameter values are p = 1.662 and Rq = 2.14, which is inside 
the endemic bubble. The period is approximately 25, which is close to the theoretically 
predicted value 2nju* ^ 24.7. 
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3.3 Response function h{y,p) = e 


-py 


This case was also studied in HU, although the bifurcation diagrams were not given in that 
paper. In this case, we can not solve in closed form the equation for the endemic equilibrium 


(3.2) 


ePy^^°’P^ = Ro{l-y{Ro,p))- 


Using (3.11), equations (3.3) and (3.4) become 


a{Ro,p) = 


-y{Ro,p) 

1 - y{Ro,p) 


and 


P{Ro,p) = -py{Ro,p). 


Implicit differentiation of (3.11) gives 


yRoiRo,p) = 


^-y{Ro,p) 


Ro{l + p{l - y{Ro,p)))' 


from which we may now compute 
<^Ro{Ro-,p) = “ 

and 


f^RoiRo^p) = 


{1 - y{Ro,p)) Ro {1 + p{l - y{Ro,p))) 
p{l-y{Ro,p)) 


Ro (1 +p(l - y{Ro,p)))' 


For r = 10, solving (3.6) becomes equivalent to solving 


-10(1 - y{Ro,p)fp^ + 11 - y{Ro,p) = 0, 


or 


P = 


1 ^10(1 “ 2 /(^ 0 ,p))(ii - y{Ro,p)) 


10 


{l-y{Ro,pW 


(3.11) 


(3.12) 
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Using this expression for p, we may solve (numerically) the equation 


a (3 cos(10\/— a^) = 0 


for y, and the result is 


y ^ 0.2703 


which, when substituted into (3.12) gives 


p* ^ 1.6617 


and then (3.11) gives 


Rl ^ 2.1474. 


We then compute 


a* -0.3704, (3* ^ -0.4491, w* 0.2540. 


Using these values and further implicit differentiations of (3.11) (we omit the details) yields 


at ^ 0.0503, at ^ -0.0190. 


Finally, the hrst Lyapunov coefficient has the value 

Kl ^ -0.4906 


so that 


£ = sgn(at/iFt) = +1, 7] ^ —0.1025 (p — p*). 


Based on Theorem 3.1, we therefore predict the following for p near p* when h{y,p) = e 
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(i) Suppose p < p*. Then varying Rq near in (1.7) will not generate Hopf bifurcation 
from the endemic equilibrium, hence there will be no endemic bubble. 


(ii) Suppose p > p*. Then varying Rq near Rq in (1.7) will generate an endemic bubble 
via Hopf bifurcation at Rq = Rq followed by a reverse Hopf bifurcation at Rq = Rq, 
with Rq and Rq near Rq, and 


Rq- Rq^ - 1 {p-p*) ^ AAm^Jp-p*. 

<^4 


Resuts of numerical simulations for this case are illustrated in Figures and IT 


4 Conclusions 

In this paper, we have performed a detailed theoretical analysis of a degenerate Hopf bifur¬ 
cation in parametrized families of RFDEs where the degeneracy arises from a violation of 
the eigenvalue crossing condition. Although our detailed computations have been performed 
for scalar RFDEs, there are no complications other than those involved with cumbersome 
notation in extending these results to systems of RFDEs. We have paid particular attention 


to the cases where the RFDE is a delay differential equation of the form (2.18), since these 


cases are quite ubiquitous in the literature, and many important phenomena in nature are 


modelled using such equations. Since we give precise conditions on the parameters of (2.18) 


to characterize the degenerate Hopf bifurcation, we expect that our paper will be a valuable 
contribution to many researchers who model phenomena using delay differential equations. 


As an application of our results, we have considered the SIS model (1.7) which was 
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studied in m- We have shown that the phenomenon of endemic bubbles which had been 
reported in m originates in a degenerate Hopf bifurcation such as the one studied herein. 

It is worth noting that the formula for the Lyapunov coefficient in the appendix includes 
the parameters a*, (3* and r. It is conceivable that by varying r, we could achieve an even 
higher-order degeneracy where both the crossing condition and the first Lyapunov condition 
are violated. This is work in progress. 


A First Lyapunov coefficient 


For the delay-differential equation (2.18), formulae (2.15), (2.16) and (2.19) yield the Lya¬ 


punov coefficient Ki in terms of the coefficients f{j,k) in (2.19) and a*, (3*, oj* and r. The 
following result was computed using the symbolic computation software package Maple: 
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Ki = 




3(1 — a*r)/(3^o) + 


3a*^T-a*^+p*‘^-3a* 

/3* 


2 , 1 ) 


2a*^T+a* l3*^T-2a*^+2a* 

- ^ - 1 Ja,2) 


I r, / 6a*^T-9a* 0* T-2a*^+2l3*^-6a*+9l3*\ r2 

(a*+/3*)(4a*-5/3*) 7 •^(2,0) 




-2 


18a*^T-33a*^0* t+ 9a* /3*^r-10a*3+7a*^/3*+10a* /3*^-7/3*^-18 q:*^+ 33a* 0*-9 0*^ 

{a*+0*){4:a*-5 0*)0* 

(a*-0*){ea*^T-9 a* 0* T-Qa*'^+a* 0*+7 0*^-9 a*+9 0*) 


2 , 0 )/(l,l) 


(a*+/3*)(4o*-5/3*)/3* 


/( 2 , 0 )/( 0 , 2 ) 


+ 


/ (o*-/3*)(4«*3t-10«*^/3* t+o* /3*2t- 40*3+2a*2/3*+3o* /3*2-3/3*3-4 q*3+10o* 0*-0*^') \ ^ 

1 0*^{a*+0*){4a*-5 0*) ) ■1(1,1) 

8 T0*3+80*^/3* T-32a*^0*^T+l9a*^0*^T-9a* 0*‘^T-8a*^-Sa*‘^0*+36a*^0*^ 

0*3(a*+0*){ia*-5 0*) 


+ 


a*^0*3_28a* 0*“^+? 0*^-8 a*‘^-8 a*^ 0*+32 a*'^ 0*^-19 a* 0*^+9 0*“^ 
/3*3(o*+/3*)(4o*-5/3*) 


1,1 


)/( 0 , 2 ) 


r, /4o*"‘T+4o*3/3*T-13o*2/3*2r+2o*/3*3r-4o*'*-4o*3/3*+15Q*3/3*2+4o*/3*3 

^ /3*2(o*+/3*)(4o*-5/3*) 


/3*2(o*+/5*)(4o* —5/3*) ' •'('1 


'( 0 , 2 ) 


It follows that in order for i^i to be well-defined, we need the non-degeneracy conditions 


a* + (3* ^ Q and 4a* — 5/3* 7 ^ 0. For a generic delay-differential equation of the form (2.18), 
this coefficient iFi will be non-zero. 
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